General Disclaimer 


One or more of the Following Statements may affect this Document 


• This document has been reproduced from the best copy furnished by the 
organizational source. It is being released in the interest of making available as 
much information as possible. 


• This document may contain data, which exceeds the sheet parameters. It was 
furnished in this condition by the organizational source and is the best copy 
available. 


• This document may contain tone-on-tone or color graphs, charts and/or pictures, 
which have been reproduced in black and white. 


• This document is paginated as submitted by the original source. 


• Portions of this document are not fully legible due to the historical nature of some 
of the material. However, it is the best reproduction available from the original 
submission. 


Produced by the NASA Center for Aerospace Information (CASI) 



N84- I J705 


(♦ 


(*ASA-TA-d4<rfd2) *U JTULAcl Ai. NEbUUNL 
CBSE H VL J bY 1.AJE05 ANT THE EFFECTIVE 
VISCOSITY OF THE EONEE flANTL*. (NASA) J 'i p 
HC AJJ/BF AO 1 CSC L 080 

C3/46 


Uacla s 
42678 


IVIASA 

Technical Memorandum 84982 


Postglacial Rebound Observed by 
Lageos and the Effective Viscosity 
of the Lower Mantle 


David Parry Rubincam 


FEBRUARY 1983 


National Aeronautics and 
Space Admir.istrction 

Goddard Space Flight Center 

Greenbelt, Maryland 20771 



POSTGLACIAL REBOUND OBSERVED BY 
LAGEOS AND THE EFFECTIVE VISCOSITY 
OF THE LOWER MANTLE 

David Parry Rubincam 

Geodynamics Branch, Code 921 
Goddard Space Flight Center 
Greenbelt, MD 20771 

ABSTRACT 

Postglacial rebound appears to have been observed gravitationally by the Lageos satellite. 
Sixty-four observations of the orbital node made over a five year time interval reveal an acceleration 
of (-8,1 ± 1 .8) X 10" 8 arcseconds day" 2 due to a source which is not presently modeled in the 
GEODYN orbit determination computer program. This acceleration cannot be explained by the 
ocean tide with 18.6 year period, assuming it to be an equilibrium tide. Instead it seems to be due 
to postglacial rebound, which changes the J 2 coefficient in the spherical harmonic expansion of 
the earth’s gravitational Field at the rate of (-8.2 ± 1 .8) X 1 0~ 19 s" 1 ; this in turn accelerates the 
node. This rate does not agree with the -32 X 10" 19 s" 1 predicted by Wu and Peltier’s (1982) L2 
model, which has upper and lower mantle effective viscosities of 10 21 and 10 22 Pa s, respectively. 

It does agree well with their LI model, which gives about -10 X 10“ 19 s” 1 . Since the effective 
viscosity is 10 21 Pa s throughout the entire mantle in the LI model, the results support the conten- 
tions that (1) the effective viscosity is near 10 21 Pa s everywhere in the mantle, and (2) this rel- 
atively low value for the effective viscosity may have permitted several degrees of polar wander due 
to glaciation during the Quaternary Ice Age. 
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POSTGLACIAL REBOUND OBSERVED BY 
LAGEOS AND THE EFFECTIVE VISCOSITY 
OF THE LOWER MANTLE 


INTRODUCTION 

During the Quaternary Ice Age the changing distribution of surface loads caused by the growth 
and decay of ice sheets deformed the earth. The deformation in turn altered the earth’s gravita- 
tional field. The remnants of these effects provide information on the earth’s rheology. In partic- 
ular, ancient shorelines and the present-day free air gravity anomalies associated with the postglacial 
rebound of Laurentide Canada and Fennoscandia (see Plates 1 and 2) provide information on the 
earth’s effective viscosity. Indeed, the relationship between deformation, gravity, and effective 
viscosity is a classical subject in geophysics (see, e.g., Cathles , 1 975 ; Peltier, 1981; Peltier et ah, 
1978; Wu and Peltier, 1 982; and references contained therein.) The value of the effective viscosity 
of the mantle has, of course, important implications for mantle convection and polar wander (e.g., 
Goldreich and Toomre , 1969; Peltier, 1980; and Lambeck, 1980). If the effective viscosity of the 
mantle is too high, then neither mantle convection nor polar wander will occur. 

The purpose of the present investigation is to infer the effective viscosity of the mantle usin^, 
satellite data. It involves looking at the rate of change of the earth’s gravitational field as deduced 
from observations of the Lageos satellite. The basic idea here is that the postglacial rebound pres- 
ently occurring changes the earth’s gravitational field as well as its geometric shape. The changing 
gravitational field in turn affects the orbits of satellites. Since the rate of rebound is controlled in 
part by the effective viscpsity of the mantle, observations of satellite orbits can in principle give 
information about the effective viscosity (e.g., Rubincam, 1979, pp. 6223-6224, and O’Keefe 
et al. , 1979; see also Paddack, 1967 and Kozai , 1966.) In particular, postglacial rebound decreases 
J 2 , the second degree, zeroth order term in the spherical harmonic expansion of the earth’s 


gravitational field; the earth is becoming less flattened. The decrease in J 2 should manifest itself 
• * 

as an acceleration £2 in the node £2 of a satellite’s orbit as the node progresses along the equator. 

Satellite laser ranging data obtained over a 5 year time interval reveals that in fact the node 

of Lageos’ orbit is undergoing an acceleration due to a source which is not presently modeled in 

the orbit determination computer program. The acceleration is presumably due to postglacial 

rebound and the ocean tide with 18.6 year period; neither of the^e are contained in the program. 

A detailed analysis of the ocean tide with 18.6 year period, given in the Appendix, indicates that 

it contributes little tc £2 over the time interval considered, assuming the tide to be an equilibrium 

one. After removing the tidal signal an acceleration still remains of £2 - (-8.1 ± 1.8) X 10" 8 

arcseconds day" 2 which is assumed to be due to postglacial rebound; this means that the rate at 

which J 2 is changing with time is J 2 »• dJ 2 /ai = (-8.2 ± 1 .8) X 10" 19 s _1 . This value is about half 

■» 

that found by Yoder et. al. (1983), who also investigate J 2 from Lageos observations. 

• 

The observed value for J 2 is compared to the values predicted by the LI and L2 earth models 
of Wu and Peltier (1982). Both of these realistic models are based on the Maxwell rheology and 
fit the Laurentide gravity anomaly and ancient shoreline data fairly well. The LI model has an 
effective viscosity for the lower mantle of 10 21 Pas (10 22 P), while the L2 model has a IQ 22 Pa s 
lower mantle. An upper limit on J 2 for LI and L2 is estimated from modeling the postglacial 
rebound of Laurentide Canada, which is where the major ice sheet of the Quaternary Ice Age rested. 
The Laurentide ice sheet is assumed to be a surface mass layer in the shape of a spherical cap whose 
mass waxes and wanes according to the ramp functions shown in Fig. 3. A lower limit on J 2 for 
the LI and L2 models is obtained from assuming that glaciation in other parts of the world, 
especially Antarctica, can reasonably be expected to give a total effect on J 2 which is 3 times 
larger than that due to the Laurentide ice sheet alone. A (weakly) preferred value for J 2 between 
the uppei and lower limits is 5/3 that of the Laurentide ice sheet. 
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These considerations lead to the following results. For model LI, -17.9 X 10" 19 s" 1 < \ 2 < 
-5.9 X 10" 19 s" 1 , with a preferred value of -9.8 X IQ -19 s’* 1 . These values agree quite well with 
the observed value of (-8.2 ± 1.8) X 10" 19 s'* 1 . For model L2, -58.4 X 10“ 19 s' 1 < J 2 <-19.2 X 
1 O’" 1 9 s" 1 , with a preferred value of -32.0 X 10" 19 s -1 . This model makes J 2 decrease much too 
fast in comparison to the observed value. Hence of the two, the LI model with its 10 21 Pa s lower 
mantle effective viscosity is preferred to the L2 model with its 10 22 Pa s lower mantle. r his 
indicates that there is little difference in the effective viscosities of the upper and lower nantle, 
a result supported by the recent studies of Yuen et. al. (1982) and Peltier and Wu (1983). 

Also, a 10 21 Pa s mantle allows a significant amount of polar wander due to glaciation - perhaps 
several degrees worth over the last few million years ( Nakiboglu and Lambeck, 1980, 1981 ; 
Sabadini and Peltier, 1981; and Sabadini et al., 1 982a, 1982b). 

LAGEOS NODAL ACCELERATION 

Lageos was launched into orbit on 4 May 1976 for the purpose of measuring crustal move- 
ments, plate motion, polar motion, and earth rotation ( Smith and Dunn , 1980; Rubincam, 1982). 
The semimajor axis of Lageos’ orbit is 1.227 X 10 7 m (about 2 earth radii) and its eccentricity 
is 0.004 (a nearly circular orbit). The orbital inclination with respect to the earth’s equator is 
109.9 degrees, while the rate at which the node of the orbit progresses along the equator is 
£2 = +0.343 degrees day -1 . 

The data consist of 64 observations of nodal residuals spread over a 5 year time interval from 
1976 to 1981 (Fig. 1). These values are what remain to be explained after running the laser range 
data through the GEODYN orbit determination computer program and after empirically determin- 
ing and removing the Kj , K ? , Pj , and S 2 tides (solid earth plus ocean) from the remaining signal 
(Peter J. Dunn, private communication, 1982). These tides have periods (from the satellites’s 
point of view) of 1051, 521, 221, and 280 days, respectively. The data points of Fig. 1 show a 
linear trend plus a slight curvature. Part of the slope can be explained by assuming that the value 
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for J 2 used in the GEODYN program needs a slight adjustment. The remaining part of the linear 
trend plus the curvature are presumably due to postglacial rebound and the 18.6 year period ocean 
tide, since neither of these are modeled in GEODYN. The 18.6 year tide is the only one which can 
contribute significantly to the curvature. It is accordingly investigated next, before proceeding to 
postglacial rebound. 

This tide is due to the precession of the lunar orbit about the ecliptic. It is probably an 
equilibrium tide because of its long period ( Proudman, 1960). Analysis of tidal records indicates 
that its amplitude is in fact close to its equilibrium value ( Currie , 1976; Rossiter, 1967). Moreover, 
Agnew and Farrell (1978) find that the amplitude of the equilibrium tide on an earth with conti- 
nents is about the same for an earth with a global ocean. Hence on the basis of these studies it 
appears that the most reasonable way to handle the data is to assume that the 18.6 year tide is an 
equilibrium one with an amplitude equal to that for the global ocean tide and subtract its effect 
from the data. In this case the perturbation in the node of Lageos’ orbit is (see the Appendix) 

A£2^ = 0.224 sin £2* arcseconds (1) 

A* 

where SI is the node of the moon’s orbit measured with respect to the ecliptic and the superscript 
“0” stands for “Ocean.” A plot of (1) (Fig. 1, solid curve) shows that A Sl^ varies almost linearly 
with time over the time interval of the data points. This is due to the fact that during this time 
ASl^ is rising from a trough to a peak and to the long period of the tide. Hence the ocean tide 
contributes mostly to the slope of the data points and only a little to the curvature. 

The signal remaining after the ocean tide contribution (1) is subtracted out is presumably due 
to postglacial rebound and the slight error in the value of J 2 . Thus its functional form should be 

Aft=£2 0 + f2 0 T+ (1/2) SI T 2 
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where T is the elapsed time since MJD (Modified Julian Date) 42904, £2 0 is a constant, £2 Q is the 
slope, and £2 is the desired acceleration of the node due to postglacial rebound. A standard least 
squares fit to the modified data yields £2 Q = 0. 10 ± 0.02 arcseconds, £2 0 = (5.7 ± 0.2) X 10" 4 
arcseconds day” 1 , and 

il - (-8.1 ± 1 .8) X 10" 8 arcseconds day” 2 
= (-5.3 ± 1.2) X 10” 23 rads” 2 (2) 

Fig. 2 shows the curvature in the data with the constant, slope, and ocean tide removed. 

It remains to find J 2 from (2). This is easily done, since the rate £2 at which the node of a 
near-earth satellite’s orbit progresses along the equator is proportional to J 2 (e.g., Kaula, 1968, 
p. 174). Hence it follows by differentiation that j ? /J 0 = Sl/il, so that 

J 2 “ (-8.2 ± 1.8) X 10" 19 s”' (3) 

after using (2) and the numerical values J 2 = 1.0826 X 10" 3 ( Stacey, 1977, p. 332) and £2 =+0.343 

« i 

degrees day" 1 . This is the value for J 2 which is assumed to be due to postglacial rebound. 

EARTH MODELS 

What must be done now is to estimate J 2 for Wu and Peltier’s (1982) earth models LI and L2 
to see whether they agree with the observed value given by (3). This will be done by first consider- 
ing the effect of postglacial rebound in Laurentide Canada on the gravitational field. This will 
provide an upper bound on J 2 (remembering that it is a negative number), since Laurentide Canada 
is only a partial (though probably the major) contributor to J 2 . Past glaciations in other regions 
of the earth, such as Fennoscandia, Siberia, and Antarctica also contribute to ) 2 . Considerations 
of these other regions, based mostly on the maximum drop in global sea level, then lead to a lower 
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limit on 3 2 • Finally, a (weakly) preferred value for J 2 between the upper and lower limits for each 
of the models will be found. 

The Laurentide gravity anomaly near Hudson Bay is quite apparent on the GEM 10B free air 
gravity anomaly map (Plates 1 and 2). The GEM 10B field is based on satellite ranging, satellite 
altimetry, and surface gravity data (L erch et ah , 1981); “GEM” stands for “Goddard Earth Model.” 
The gravity anomalies were computed from the standard equation (e.g., Rapp, 1975, p. 198) 

Ag = 7 £ (fi-D £ (Cg m cos mX + S fim sin mX) P Km (cos 0) 
i B 2 m n 0 

(written here in normalized form) using the GEM 10B gravitational field coefficients (Cg m , Sg m ) 
given by Lerch et al. (1981) up to and including degree and order 36, except for the £ = 2, 4, 
m = 0 terms. For these terms the hydrostatic values C 2 q = -480.5 16 X IQ" 6 and - +1.212 X 
10“ 6 derived from Nakiboglu (1979, p. 645) were subtracted from the corresponding GEM 10B 
values. Hence the anomalies shown in Plates 1 and 2 are basically referred to the earth’s hydro- 
static figure. The underlying maps shown in the plates are the tectonic activity maps of Lowman 
(1981, 1982). 

Note that the gravity low in Laurentide Canada reaches a minimum value of about -50 X 10" 5 
m s" 2 (-50 mgal) when referred to the hydrostatic flattening. This is in contrast to the usual 
-40 X 10" 5 m s" 2 when referred to the reference flattening f - 1/298.255 (e.g., Wu and Peltier, 
1982; Peltier and Wu, 1982; Rapp , 1975, p. 210). 

The ancient Laurentide ice sheet is modeled here as a surface mass layer a ICE of constant 
density and in the shape of a spherical cap. The angular radius of the cap is a and its center is 
located at colatitude (3 . The mass M of the ice sheet varies with time t according to the equation 
M = M 0 f(t), where M 0 is a constant and equal to the maximum mass of the ice sheet and f(t) is 
shown in Fig. 3. The function f(t) represents a simplified version of the glaciation-deglaciation 
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history of Laurentide Canada (Sabadini and Peltier, 1981). The ramp functions of Fig. 3 are 


characterized by the constants Tj and T 2 . The accretion time Tj is ~ 100 000 years and the 
disintegration time T 2 is ~ 10 000 years. The time T 3 s 5000 years is the elapsed time since the 
disappearance of the last ice sheet. Also, the cycles extend about 2 million years into the past 
(e.g., Wu and Peltier, 1982, p. 480.) 


The water composing the ice sheet comes from the oceans, assuming a closed hydrologic cycle 
Hence the oceans must also be represented as a surface mass layer o oc which varies with time. 
Sabadini and Peltier (1981, p. 558) give tiie equations for a ICE and a oc , Adding them together 
to obtain the total mass layer a - o ICE + a oc yields 


o- £ o 2 (0', t) 


where 


o e (0',t) 


M 0 f(t) 
2 tt R| 


(2£-H) dFgpCcos oQ 
fi(g+l) 9 (cos a) 


^ 0 (cos O'). 


(4) 


Here P^qCcos 6') is the associated Legendre polynomial of degree 2 and order 0, O' is the angle 


measured from the center of the cap, and R E is the radius of the earth. 


The effect of this ice sheet on the present-day exterior gravitational field must now be found. 
It will be given by 

v i = (^)/ ^(0",t / )G fi (\i/,t-t')dA , 'dt' (5) 

4.00 

where Vg is the term of degree 2 in the spherical harmonic expansion of the gravitational potential 
about the axis of the cap, r is the radial distance, dA" is the element of area, and ug(i//,t-t') is the 
Green function ( Wu and Peltier, 1982, pp. 468-469) 

Gj(iM-t) = Gf (iM-t’) + Gj' (list-t') 
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where 


Gf ~ kg fi(t-t') Pg 0 (cos M 

k e 


( 6 ) 


and 


GifaM-f )-.£.£ r? e^ (, '° 


jg y/,i— l ) - ^ jy e ' PgpCcos ^ 


(7) 


Here G is the universal constant of gravitation, Gjf is the clrrtic part of the Green function, kg is 
the elastic Love number, S(t-t') is the Dirac delta function, Gg is the viscous part of the Green 
function, ^ is the angle from the load point to the point of observation, K is the number of 
relaxation modes, and the i| and sj 2 are numbers which characterize the particular earth model 
and are discussed below. The associated Legendre polynomial Pg 0 (cos i/O can be written 

X (2-5 0m )(fi-m)! 


Pg 0 (cos 0) = 2 


i*»V 

nro 


(Ki-m)! 


• Pg m (cos Q ") P Crn (cos O') • [cos mX" cos mX' 


Cm 


+ sin mX" sin mX'] (8) 


using the addition theorem for spherical harmonics (e.g., Kaula, 1968, p. 67). Here 0" and O' are 
colatitudes and X” and X' are east longitudes measured with respect to the axis of the cap. 


The elastic part (5) of the Green function is of no interest here and will not be considered 
further. Also, it proves to be slightly more convenient to use nj = 1/rj 2 and rj = l/sj rather than 
i| and s^ ; the numbers which characterize the earth are now measured in years rather than inverse 
years. Further, only the £ = 2 parts of (4), (5), (7), and (8) need be found since this investigation 
is concerned only with the second degree term in the potential. Hence setting £ - 2 in these equa- 
tions, using (8) in (7), (4) and (7) in (5), and carrying out the integration over area gives 

Pg 0 (cos O') 

where 3P 20 (cos a)/d(cos a) = 3 cos a has been used. 
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Applying the addition theorem (8) once again to express the potential in terms of spherical 
harmonics about the rotation axis of the earth rather than the axis of the cap and evaluating the 
time integral using the function shown in Fig. 3 yields 

. GM„ , /R e \ 3 

v *“ii7 J2 \~r) Pj ° (cos(,) 

where 


, = ”5- cos “(t cos2 ' 5 -t) 


1 



1 

>> 

1 



JL+ £ 


-Ndj+TjJ/Tj 


'JL+JL 

T, T 2 


j_ e ~ T V 7 f + e Wi* t 2> + Ti W) _ e -<N-l)(T 1+ T 2 )/r j 2 


/ t i + t 2 )/T _ 


-yi 


(9) 


Here V 2 is the desired second degree coefficient in the potential, 0 is colatitude measured from 
the North Pole, M E is the mass of the earth, N is the number of glacial cycles, and the superscript 
“L” on V 2 and stands for “Laurentide.” Use has been made of the expression S = (l — X N )/(1 — X), 
where S = 1 + X + X 2 + . . . + X N “ 1 , in evaluating the time integral. Taking the derivative of (9) 
with respect to time T 3 finally gives 


'" = -£ cos “(f oos 2 ^-t) 


K 

S 

j=] 

i 

\i] 


r ±+e -N wf / 

wJ 


_ T 2 Tj [ 



l-e~ Tl/T i + e l " N(T l +T 2 )+T l l/r J _ e -<N-D(T 1+ T 2 )/r 2 

» e (T l +T X „ , 



( 10 ) 


as the rate of change of J 2 due only to the Laurentide ice sheet. 


Ael 


The question now aris ;s as to wIsh values to use for the parameters appearing in (10). The 
values for M 0 , /3, and a are taken from Sabadini et al. (1982a, p. 2897) and are given in Table 1. 

The values for Tj , T 2 , and T 3 , also given in Table 1 , are estimated from Fig. 2 of Andrews and Berry 
(1978), which shows the last glaciation-deglaciation cycle of the Quaternary Ice Age, The times 
Tj ~ 1 12 000 and T 2 « 13 000 years are somewhat longer than those used by Sabadini and Peltier 
(1981), whose values of 90 000 and 10 000 years are more in accord with the period of the 
eccentricity cycle of the earth’s orbit which apparently drives the glaciation ( Hays et al. 1976). 

The longer times are used here since the last accretion-disintegration cycle is the most important; 
the earth will not “remember” much of the previous cycles for the relatively low effective viscos- 
ities considered here (Sabadini and Peltier, 1981, p. 568). The implications of using the shorter 
astronomical times will be discussed below. The Ice Age extends about 2 million years into the 
past, so N s 20 (e.g., Wu and Peltier, 1982, p. 480). The values for M E and R E come from Stacey 
(1977, pp. 331-332). 

The values for r 2 = 1/ts? and r 2 /« j = r 2 /s 2 for the LI and L2 models are derived from 
Wu and Peltier (1982, pp. 465-466) and are shown in Table 2. Both models have the same density 
and elastic properties as model 1066B of Gilbert and Dziewonski (1975). Model LI has a 120 km 
thick lithosphere, a mantle with effective viscosity of 10 21 Pa s, and an inviscid core. Model L2 
is the same as LI except that the effective viscosity of the mantle below the density discontinuity 
at 671 km depth is 10 22 Pa s. The multiple relaxation modes are due to the discontinuities present 
in the models: M0 is the fundamental mantle mode, Ml and M2 are due to the density discon- 
tinuities at 671 km and 420 km depth, respectively, and L0 and CO are due to the presence of the 
lithosphere and the core, respectively. Both the LI and L2 models fit the Laurentide gravity 
anomaly referred to the nonhydrostatic figure and bracket the relative sea level data, leading 
Wu and Peltier ( 1982) and Peltier and Wu (1982) to conclude that the effective viscosity of the 
lower mantle is between 10 21 and 10 22 Pa s. 
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Using the numerical values given in Tables 1 and 2 in (10) yield the following results, For 
model Li 

J 2 s -5.9 X 10“ 19 s“ l (11) 

wliile for model L2 

=-19.2 X 10" 19 s" 1 (12) 

for postglacial rebound due to Laurentide Canada ulone. Note that model L2 already predicts a 
rate of decrease in J 2 which is over twice the observed value given in (3), while model LI predicts 
a rate which is slightly less than the observed value. 

The rates given by (1 1) and (12) must be considered upper bounds on J 2 , since Laurentide 
Canada is not th* only region to have had an ice sheet. Other parts of the globe underwent giaeia* 
tion n well o .u will also contribute to J 2 , A rough estimate of their contributions will be considered 
next. 

The maximum drop in global sea level was probably about 80 m and may have been as high as 
120 m (e.g., Andrews and Berry, 1978, p. 210). It can easily be shown that the Laurentide ice sheet 

accounts for only 48 m or so of this drop, assuming that the ice sheet had the mass given in Table 1 ; 

* 

other ice sheets must account for the remaining 32 to 72 m. Certainly a lower limit on J 2 can be 
obtained by assuming that 72 m was due to an ice sheet comparable to the Laurentide ice sheet 
situated in Antarctica over the South Pole. (There are in fact large Antarctic anomalies which may 
be at least partially explained by glaciation; see Plate 2). This maximizes the effect on J 2 not only 
by maximizing the allowable mass but also by placing it at the pole (the angle /3 in (10) is now 180 
degrees rather than 25 degrees). Use of (10) gives this hypothetical ice sheet an effect on J 2 which 
is about twice that of the Laurentide ice sheet; thus the total effect on J 2 is about 3.T£ so that 

-17.9 X 10" 19 s~ l < J 2 (13) 


il 


for model LI and 


-58.4 X 1CT 19 s" 1 < J 2 (14) 

for model L2. These calculations are rough since they assume an Antarctic ice sheet with the same 
radius as the Laurentide ice sheet (but with 1.5 times the mass) and the same history; they also 
assume that the effect of each ice sheet may be calculated independently of the other, i.e. there 
is no coupling between the two. However, only rough estimates are desired because the uncertain- 
ties of the sea level data make further refinements not worthwhile. 

Preferred values for J 2 for models LI and L2 will now be found. The Fennoscandian, 
Siberian, and other ice sheets (excluding Antarctica) had a combined mass which was about 42 per 
cent that of the Laurentide ice sheet ( O’Connell, 1971, p. 308) and accounted for about 20 m of 
the global drop in sea level. Hence by adding in the 48 m or so due to the Laurentide ice sheet 
about 48 + 20 = 68 m of the probable 80 m drop in global sea level can be explained. Thus by 
assuming that these other ice sheets in fact account for 32 m of an assumed 80 m total drop and 
by lumping them together into one ice sheet with the same radius and colatitude as the Laurentide 
ice sheet but with 2/3 the mass yields J 2 = (5/3) j^. So for model LI 

J 2 = -9.8 X 10" 19 s” 1 (15) 

and for model L2 

J 2 =-32.0 X J O -1 9 s" 1 . (16) 

These values are only weakly preferred due to the limitations of the calculations as discussed 
above and because a possibly substantial amount of deglaciation in Antarctica has been ignored 
(Sabadini et al., 1982a, pp. 2897-2898). Even so, the value (15) agrees well with the observed 


value (3), while (16) is a factor of 4 too large. Hence of the two models, LI is preferred to L2. 


DISCUSSION 

It has been assumed here that the curvature evident in the nodal residuals is due mainly to 
postglacial rebound. There may of course be alternative ways of explaining the data. One such 
alternative which immediately comes to mind is the 18.6 year ocean dde; perhaps it is not an 
equilibrium tide as assumed here. In fact, Sanchez (1979) speculated that this tide may have an 
amplitude several times that of its equilibrium value in order to explain a polar motion component 
with the same period found by Markowitz (1979) from polar motion data. 

The data shown in Fig. 1 were accordingly reanalyzed to investigate this possibility by assum- 
ing that the functional form of the nodal residuals is now given by 

A£2 = £2 0 + £2 0 T + I &T 2 + Aj sin . 

M 

Here the last term represents an 18.6 year ocean tide which is unconstrained in amplitude but 

constrained to having the equilibrium phase. In this case solving for £2 0 , £2 0 , £2, and Aj gives 

Aj = 0.42 ± 0.15 arcseconds and £2 = (-5.9 ± 2.5) X 10" 8 arcseconds day" 2 , which would indicate 

that the tide has about twice its equilibrium amplitude given by (1). But this still does not explain 

the quadratic behavior of the residuals; the new value for £2 is not significantly changed from (2) 

and still agrees fairly well with Wit and Peltier’s (1982) LI model. Of course the ocean tide could 

A ^ 

be unconstrained in both amplitude and phase by adding a term, of the form A_ 2 cos £2 to the 
above equation and solving for all of the coefficients, including A 2 . In this case £2 = (4 ± 100) X 10" 8 
arcseconds day" 2 , A 1 = 0.44 ± 0.23 arcseconds and A 2 = -0. 12 ± 1 .3 arcseconds, so that £2 is not 
significantly different from zero and no relaxation of the earth has been observed. However, it is 
unreasonable to solve for the amplitude and phase of a tide with an 18.6 year period with only 5 
years worth of data. In any case the ocean tide would have to depart drastically from its equilib- 
rium value in order to explain the Lageos observations. 
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A different ocean tide might appear to be another way of accounting for the T 2 behavior of 
the points in Fig. 1 ; but (A3) jidicates that there is no tide of suitable period and amplitude to 
do so. 

The deceleration of the earth’s spin by tidal friction does not account for it either, as may be 
seen from the following considerations. For a homogeneous, liquid earth the hydrostatic flattening 
factor f H is proportional to co 2 to a first approximation, where co is the angular speed of the earth 
(e.g., Lyttleton, 19S3, p. 38). Further, elementary considerations show that also to a first approx- 
imation J 2 is proportional to f H (e.g., Kaula , 1968, pp. 68-69). Hence by differentiation i 2 /J 2 ~ 
2co/co, where cb is the acceleration of the rotation speed due to tidal friction. Using the values 
co — —6 X 10" 22 rad s” 2 (Stacey, 1977, p. 99) and co = 7.29 X 10" 5 rad s" 1 ( Stacey, 1977, p. 332) 
in this equation give J 2 == -2 X 1G“ 20 s" 1 , which is a factor of 40 smaller than the observed rate 
given by (3). Thus even if the earth were completely fluid, the observed relaxation of the earth 
could not be due to the slowdown of its rotation rate by tidal friction. 

There may be some other alternative, as yet unthought of, which will explain the data; but 
thus far postglacial rebound does quite well. Moreover, postglacial rebound predicts the continued 
quadratic behavior of the residuals in the future; its parabolic nature should become more and more 
evident as data is taken future years. Hence the relaxation of the earth is an hypothesis which 
is easily tested. 

The data show good agreement with Wu and Peltier’s (1982) LI model for the parameter 
values adopted here. It has been mentioned that the values of the accretion time Tj = 112 000 
years and disintegration time T 2 = 1 3 000 years are longer than the 90 000 years and 1 0 000 years 
used by Sabadini et al . (1982a, 1982b) and Sabadini and Peltier (1981). Use of these shorter times 
give a preferred J 2 = -12.5 X 10 -19 s" 1 for the LI model and -36.7 X 10" 19 s -1 for the L2 model, 
as compared to the observed (-8.2 ± 1.8) X 10“ 19 s -1 . Thus the shorter times, which appear to be 
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more in accord with the deep-sea sediment core data (e.g., Hays et al., 1P76; Sabadini and Peltier, 


1981, p. 566), worsen the agreement slightly, but not seriously, between the LI and observed 
values. 

The results given here can be compared to other recent studies. Yoder et al. (1983) have 

t 

independently analyzed observations of Lageos’ node and have concluded that J 2 is about 
(-16 ± 3) X 10 -19 s" 1 . This figure is twice as large as that found here, and falls between the 
preferred values (15) and (16) for the LI and L2 models, arguing for a lower mantle effective 
viscosity between 10 21 Pa s and 10 22 Pa s. Also, Yuen et al. (1982) inferred the effective viscosity 
of the lower mantle from the observed nontidal acceleration of the earth’s rotation. They find the 
effective viscosity to be between 1 and 4 X 10 21 Pa s, in good agreement with the results given here. 
Finally, Peltier and Wu (1983) also calculate from rotational acceleration that it should be between 
1.6 and 6 / 10 21 Pa s. On the whole it appears that the effective viscosity of the lower mantle is 
indeed between 10 21 Pa s and 10 22 Pa s, with the present investigation supporting the lower value. 

It is interesting to note that the L2 model has J 2 decay more rapidly than the LI model, even 
though its effective mantle viscosity is a factor of 10 higher than that of the LI model. This seem- 
ingly paradoxical behavior is due to the discontinuities present in these realistic models ( Wu and 
Peltier , 1982, p. 475-476). This behavior also shows that previous attempts by Rubincam (1979) 
and O’Keefe et al. (1979) to derive the effective viscosity of the lower mantle from satellite data 
using an incompressible, homogeneous, viscous earth in the manner of Darwin (1879) are too sim- 
plistic. It can be shown that the classical Darwinian model applied to the Lageos data give an effec- 
tive viscosity of about 10 23 Pa s — far from the more reasonable value of 10 21 Pa s. Hence the data 
cannot be inverted to give a unique effective viscosity. Instead realistic models must be constructed 
and then tested against observation. 



It is also of some interest to compare the glacial flattening to the total nonhydrostatic tlatten- 
ing of the earth, Wu and Peltier (1982) do not give quite enough information in their Table 8 to 
compuv? for the LI model, but do in their Table 9 for the L2 model, which gives = 1.1 X 10" 6 
by (10). (Both the LI and L2 models agree with the gravity anomaly data referred the reference 
flattening and should give about the same value.) The preferred value is J^” = 1.8 X 10" 6 with a 
maximum possible value of 3.3 X 1 0 -6 . The total nonhydrostatic value AJ 2 for the earth is about 
8. 174 X IQ' 6 , as may be derived from Table 3 of Nakiboglu (1979). So the squashing of the earth 
by the glaciers of the Quaternary Ice Age accounts for between 13 per cent and 40 per cent of 
the present-day nonhydrostatic part of J 2 , with a preferred value of 22 per cent. Wang (1966) had 
suggested that glaciation might account for all of AJ 2 , but McKenzie (1966), Kaula (1967), and 
O’Connell (1971) indicated that glaciation accounts for an amount closer to the percentages given 
here. (For further discussion of AJ 2 , see Goldreich and Toomre , 1969). 

The support for the LI model indicates that there is little, if any, jump in the effective viscosity 
across the density discontinuity at 671 km depth ( Wu and Peltier , 1982; Peltier and Wu, 1982). Also, 
an effective viscosity of 10 21 Pa s for the whole mantle has important implications for polar wander. 
Sabadini et al. (1982a, 1982b), Sabadini and Peltier (1981), and Nakiboglu and Lambeck (1980, 
1981) find that glaciation on an earth with this value for the effective viscosity of the mantle could 
have caused several degrees of polar wander over the last few million years, and can explain the 
current movement of the pole towards Canada. This supposes, of course, that the linear Maxwell 
rheology assumed here applies to tectonic stresses and time scales, which may not be the case 
(e.g., Weertman , 1978; but see Wu and Peltier, 1982, p. 476). McAdoo (1982), for instance, finds 
that non-Newtonian flow best explains the geoid over subducting slabs. 
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SUMMARY 


The node of Lageos’ orbit is undergoing an acceleration not presently modeled in the orbit 
determination computer program. It appears to be caused by the postglacial rebound of the eaith. 
An effective viscosity of IO 21 Pa s for the whole mantle is consistent with the observed acceleration. 
This value for the effective viscosity allows rapid polar wander due to glaciation. 
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APPENDIX 


The problem addressed here is to find the effect of the 18.6 year global equilibrium ocean tide 
on the node S2 of Lageos’ orbit as given by (1). 


The amplitude of the global equilibrium ocean tide is (Agnew and Farrell, 1978, equation 2.5) 


$2 


1 +k 2~ h 2 

1-^ (l+k^-h') 

5p E 


g 


where U 2 is the tide-raising potential at the earth’s surface, p w = 1022 kg m“ 3 is the density of sea 
water, and k 2 , h 2 , k 2 , h 2 are the usual Love numbers. The ocean tide will produce a disturbing 
potential (Hendershott, 1972, equation 5) 


(l+k 2 )ga 2 f 2 



while the solid earth tidal disturbing potential is 


k 2 u 2 



where a 2 = 3p w /5p R . Hence the ratio between the two is 


3 (l+k 2 ) (l+k 2 -h 2 jp w 


5k» 


’Pw' 

, Pe , 


l--f-^ (l+k 2 ~h 2 ) 


= 0.22 


Pe 


(Al) 


for the values k 2 = 0.30, k 2 = -0.31, h 2 = 0.61, and h 2 = -1.00 ( Lambeck , 1980, pp. 13-14). 
Thus the ocean tide effect is 22 percent that of the solid earth tide. Hence once the effect of the 
solid earth tide on the node of Lageos’ orbit is found, the ocean tide effect immediately follows. 

19 


PRECEDING PAGE BLANK NOT FILMED 


QRIQAftiAL PAGE jg 
OF POOR QUALITY 


The solid earth tide is found next. Goad ( 1 977) argues that for the sake of easy integration 
the solid earth tidal disturbing potential should be expressed in terms of the satellite’s Keplerian 
elements (a, e, I, co, £2, M) referred to the earth’s equator and the moon’s Keplerian elements 
$*, A*, £1*, M*) referred to the ecliptic. (See also Musen and Estes. 1972.) This will be 

done here. However, Goad ’s (1977, equation 1.7) expression will not be used here since it contains 
some errors, as Goad himself realized; instead a slightly different derivation is sketched below. 


Kaula (1964) shows how to obtain the solid earth tidal disturbing potential U SE expressed in 
terms of the satellites’s and moon’s equatorial Keplerian elements from the tide-raising potentir 1 



y y (2-6 

U (*-m )! 1 


m=o 


om 


>wa)Y gmi( 0*,x*) 


(A2) 


expressed in terms of the spherical harmonics Yg ml (0,X) = Pg m (cos 6 ) cos mX, Y gm2 (0,X) = 

Pg (cos 0 ) sin mX which refer to the earth’s equatorial system. The lunar equatorial spherical 

)(<>{< A 

harmonics Yg mi (0 ,X ) may be expressed in terms of the ecliptic colatitude 6 and longitude X 

using (Messiah, 1963, pp. 1073-1074) 



where the coefficients Ag mist depend on the elements of the rotation matrix between the equa- 
torial and ecliptic coordinate systems. After substituting this expression in (A2) and thereafter 
following Kaula’s (1964) treatment gives U SE = X] Ug^pq^j, where 
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J - cos 

_ v Cmhj + v Cspq “ Y 

+ cos 

. v 2mhj ~ v Cspq + 2 ] \ 

1 -sin 

. v Cmhj + v fispq ~ Y 

- sin 

. v 2mhj ~ v Cspq + *2 J / 

J -sin 

. v 2mhj + v fispq “ Y 

+ sin 

^ v £mhj “ v £spq + Y ] 1 

\ + cos 

_ v 2mhj + ''fispq-^+Of^ 

+ cos 

__ v 2mhj ~ v fispq + ^"i) Y ] / 


£-m even, fi-s even 
J2~m even, S~s odd 

J2-m odd, £-s even 
£-m odd, 2-s odd 


where 


(A3) 


’.mhj = (^~2h) co + (£-2h+j) M + m£2, 

A>je Aa, A |i< A* 

Vg spq = (£-2p) co + (g-2p+q) M* + s£2 . 

Here k g is the solid earth Love number of degree £, M m the mass of the moon, R E the equatorial 
radius of the earth, and G the universal constant of gravitation, 


For the 18.6 year tide 2-2, m = 0, s = 1 , p = 


ySE _ GM m 
u 2011010 


- — k, 

R„ 2 


Rr 
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' r e 

TSf 
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1 , q - 0, h = 1 , and j = 0, so that (A3) becomes 
R 201 R 211 ^ ^ 210 ^ e -^ ^ 210 ^ e>f ^^ 


{ ^ 20111 3 * n ^ “ A 20 j 12 COS 12 | 


(A4) 


The coefficients A 20ill and A 201i2 can easily be found from the addition theorem for spherical 
harmonics (e.g., Kaula, 1964, p. 67), since m = 0 in the above equation. They are A 20111 = 0 and 
A 201 12 “ s * n e cos e > w h ere e is the obliquity of the ecliptic. Using these values in (A4) and noting 
that G 210 (e) = U 210 (^*) ~ ^ because of the nearly circular orbits of the moon and Lageos give 


U SE 

U 2011010 


GM 


M 


Rr 




whi- h varies sinusoidally with time with a period of 1 8.6 years. Substituting this expression in 
the Lagrange planetary equation for df2/dt (e.g., Kaula , 1968, p. 166) and integrating with respect 
to time yields 




AS2 £* 


-9 GM m k 2 


4VGM E a R e 



» . 'O'# A* . . A# 

cos I sin I cos I sin e cos e sin £2 


for the effect of the 1 8.6 year solid earth tide on the node of Lageos’ orbit. Multiplying this by 

A* 

0 22 as dictated by (Al) and using the numerical values e = 23.44 degrees, I - 5.15 degrees, etc. 
(e.g., Kaula, 1968, p. 186) yields tho effect of the 18.6 year global equilibrium ocean tide as 
given by (1). 
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Plate 1. The GEM 10B free air gravity anomalies referred to the earth’s hydrostatic figure. 

Van der Grinten projection. 




Polar orthographic projection. 













FIGURE CAPTIONS 


Plate 1. The GEM 10B free air gravity anomalies referred to the earth’s hydrostatic figure. 

Van der Grinten projection. 

Plate 2, The GEM 10B free air gravity anomalies referred to the earth’s hydrostatic figure. Polar 
orthographic projection, 

Figure 1 . Plot of the residuals in the node of Lageos’ orbit after modeling the orbit with the 
GEODYN program; also the 18.6 year global equilibrium ocean tide (solid curve). MJD ® Modified 
Julian Date. 

Figure 2. Plot of the residuals in the node after removing the 18.6 year ocean tide perturbation, 
the constant, and the slope, leaving only the curvature presumably due to postglacial rebound. 

Figure 3. Assumed ice sheet history for Laurentide Canada. Tj is the characteristic accretion time 
and T 2 is the characteristic disintegration time, while T 3 is the elapsed time from the end of disin- 
tegration to the present. 
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Table 1 

Preferred Values for the Laurentlde Ice Sheet Parameters 


Parameter 

Symbol 

Value 

Radius of Ice Sheet 

a 

15 Degrees 

Colatitude of Center 
of Ice Sheet 

fi 

25 Degrees 

Maximum Mass of 
Ice Sheet 

M 0 

1.8 X 10 19 kg 

Characteristic 
Accretion Time 

Ti 

112 000 Years 

Characteristic 
Disintegration Time 

t 2 

13 000 Years 

Elapsed Time Since 
the End of Disintegration 
to the Present 

*r» 

*3 

5000 Years 


34 












TABLE CAPTIONS 


Table 1, Preferred values for the Laurentlde Ice sheet parameters, 

Table 2, Valuos for K 2 * 1/r 2 and r 2 « 1/s 2 for models LI and L2. Derived from Tables 8 and 9 
ofW u and Peltier (1982), Only those values which contribute significantly to J 2 are shown here, 


